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Abstract 


In the field of Turbomachinery, apart from design point of view, viscous flow pre- 
dictions carry greater importance for the analysis of off-design conditions: which 
would be possible by computing unsteady compressible Navier-Stokes equation. In 
the present work, viscous flow past cascades have been studied. 

An unsteady two-dimensional viscous compressible Navier-Stokes equation 
is being solved in this context by Finite Yolume(FV) principle using cell centered ap- 
proach. A higher order upwinding scheme (MUSCL interpolation) along with Steger 
and Warming Flux Vector Splitting has been implemented. A time-step marching is 
carried out with Rational Runge-Kutta(RRK) time integration method. Elliptically 
generated grid is used to compute both Euler and Navier-Stokes equation. 

The flow has been analysed for NACA0012 as well as DCA cascade for 
various different inflow Mach numbers in subsonic and transonic condition. Finally, 
various types of finite volume schemes have been theoretically analysed by spectral 
method and their efficacy have been discussed. 



Chapter 1 


Introduction 


1.1 Role of CFD in Turbomachinery 

The Physical aspects of any fluid flow are governed by three fundamental principles 
based on conservation of mass, momentum and energy. These fundamental prin- 
ciples can be expressed in terms of mathematical equations, which in their most 
general form are usually partial differential equations. Computational Fluid Dy- 
namics (CFD) is the science of determining a numerical solution to these governing 
equations of fluid flow whilst advancing the solution through space or time to obtain 
a numerical description of the complete flow field of interest. 

The governing equations in fluid dynamics, the unsteady Navier-Stokes 
equations, have been known for over a century and half. However, the analytical 
investigation of reduced forms of these equations is still an active area of research 
as is the problem of turbulence closure for the averaged form of the equations. 

Experimental fluid dynamics plays an important role in validating and de- 
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lineal, mg the limits of the various approximations to the governing equations. The 
wind tunnel, for example, as a piece of experimental equipment, provides an effective 
means of simulating real flows. Traditionally this has provided a cost effective alter- 
native to full scale measurement. However, in the design of equipment that depends 
critically on the flow behaviour, for example the aerodynamic design of an aircraft, 
full scale measurement as part of the design process is not practical and sometimes 
not feasible. This situation has led to an increasing interest in the development of 
a numerical wind tunnel. 

The role of CFD in engineering predictions has become so strong that today 
it may be viewed as a new third dimension of fluid dynamics, the other two dimen- 
sions being the above stated classical cases of pure experiment and pure theory. 

The development of more powerful computers have furthered the advances 
being made in the field of computational fluid dynamics. Consequently CFD is now 
the preferred means of testing alternative designs in many engineering companies 
before final experimental testing takes place. 

1.2 Aim of the Thesis 

A computational study utilizing higher order upwinding scheme is being proposed 
to solve two dimensional unsteady viscous compressible Navier-Stokes equation for 
flow past cascade. The aim is to study the effect of gust on the performance of the 
cascade. Here the gust is unsteadiness of the oncoming flow which is due to the shed 
wake of the previous stage. 

Thus the computational study of such a complex flow phenomena would 
require solving the three dimensional unsteady Navier-Stokes equation. Here in this 
project, we are going to study the two dimensional unsteady compressible flow past 
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cascade. The understanding and prediction of such viscous flows are of great interest 
and importance in turbomachinery design and analysis. 

Here for the Xavier- Stokes solver we have used the flux vector splitting 
method utilising higher order upwinding in a finite volume formulation. The higher 
order upwinding is based on MUSCL interpolation and can be used for analysis of 
aeroacoustic problem in cascade as well. 

1.3 Thesis Overview 

For the basics of unsteady compressible flows and fundamentals of CFD, we have 
followed Fletcher [1], Hirsch [9], Anderson [15], Courant &: Friedrichs [16], White [17] 
and David [18]. 

The development of theoretical methods to predict unsteady flows in tur- 
bomachinery is a difficult task and determining the time-dependent two or three 
dimensional flows of a viscous compressible fluid through a geometric configuration 
is of enormous complexity. An unsteady viscous phenomenon in case of a tur- 
bomachinery include intricate shock patterns, viscous separation and many other 
physical problems. Because of the complexity of the unsteady phenomena occurring 
within the machine, design and development approaches have been based generally 
on steady-state theoretical aerodynamics along with the empirical correlations to 
account for unsteady effects. 

For unsteady viscous compressible flows, explicit schemes are more accurate 
in predicting the flow solution as that of implicit schemes. For numerical stability, 
however, the time-integration step size used in integration must be fairly small, 
while allowable time step is governed by CFL number. In this thesis we have used 
these aspects to solve unsteady viscous two dimensional compressible Navier-Stokes 
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equation for flow past cascade. 

There are three major components of this thesis problem. First and foremost 
is the grid generation, which is needed to study flow past a complex body. Here, grid 
generation has been carried out by two methods: First is an Algebraic and second is 
an Elliptic grid generation method. The topic grid generation is very well described 
by Fletcher [1]. Thompson et al [2] give details about an elliptic grid generation 
technique for internal as well as external flow domain and this has been followed in 
this thesis. 

The second and the most important component is the Navier-Stokes solver. 
The formulation and flux vector splitting ideas are well described in Shu et al [3] and 
Hirsch [9]. The flux vector splitting scheme given by Steger and Warming [4] gives 
better efficiency in calculations of eigenvalues and eigenvectors while constructing 
the split fluxes. This flux vector splitting has been implemented in the solver. Mul- 
tidimensional splittings given by Deconinck [19] give the mathamatical approach 
of the splitting schemes. The MUSCL interpolation scheme during discretization 
has been used for better accuracy, which is also given in Hirsch [9]. For time-step 
marching, the Rational Runge Kutta method by Morinishi and Satofuka [10] is being 
used. The boundary conditions are implemented as given in Ravichandran [11] and 
Mandal [7], Finally, the viscous terms calculation and full Navier-Stokes solution 
is carried out following Shu et al [3] and Kallinderis & Baron [13]. For the com- 
putations, geometries are given in Fottener [6]. For NACA0012 straight cascade, 
the pressure variation along top and bottom surfaces and Mach number distribution 
has been compared with Mandal [7]. While, for the DC A cascad, the experimental 
results given in Fottener [6] for the pressure coefficients on surface, Mach number 
variation and the velocity profile close to aerofoils have been compared with com- 
puted solutions. 
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Also, we have performed the spectral analysis of flux vector splitting finite 
volume method for compressible flows in general. Sengupta and Sengupta [14] gives 
the details about spectral analysis for an incompressible flows which has been ex- 
tended here. We have analysed various flux vector splitting schemes, for their phase 
errors and numerical dissipation, which are given in Hirsch [9]. 



Chapter 2 


Governing Equations for 
Compressible Flows 


2.1 Introduction 

Most compressible flows of practical importance that require consideration of the 
full Navier-Stokes equation are turbulent. The external compressible flows at large 
Reynolds numbers, viscous and turbulence effects are only significant close to solid 
surface in the absence of massive separation. Consequently, for such flows occasion- 
ally dissipative terms associated with the direction normal to the surface are as Thin 
Layer Navier-Stokes (TLNS). 

However, for the transonic flow past cascades, the complicated shock bound- 
ary layer interaction produce local region of separated flows and it is appropriate 
to use full Navier-Stokes equation to capture the phenomenon accurately. If the 
flow is locally supersonic, embedded shock waves are likely to be present. If these 
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weak, accurate solutions can be usually obtained without modification to the basic 
computational algorithm, other than the inclusion of additional numerical dissipa- 
tion. 

This thesis discusses a Finite Volume Method(FVM) for discretisation of 
full Navier-Stokes equation. Finite Volume Method(FVM) is called as local method 
since the discretised equations are sparse. The finite volume method has additional 
advantage of discretising directly the conservative form of governing equation. This 
implies that the discretised equation preserves the conservation laws, reducing nu- 
merical errors. 

The Flux Vector Splitting(FVS) scheme, that has been used here, can be 
applied to the conservation law form which allows the shock waves to be captured as 
weak solutions to the governing equation and circumvents the difficulty of applying 
shock fitting techniques to arbitrary flows. 

Since the fluxes and geometric terms are evaluated only at the interfaces, 
the higher order MUSCL-type interpolation in space plays a major role in solution 
accuracy. Higher order Runge-Kutta methods for time integration similarly plays 
an important role in the selection of proper time step. A two-stage Rational Runge- 
Kutta(RRK) method has been implemented for the flow past cascade. 

The application of boundary conditions is also crucial to have a correct 
solution, especially, for the internal transonic flows. We have used characteristics 
based boundary conditions. 
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2.2 Finite Volume Formulation for Compressible 
Flow Equations 

2.2.1 Introduction 

The computation of inviscid or viscous compressible flow past cascade is difficult 
because of the presence of shock with the strong gradients of finite thickness apart 
from the shear layer close to the body. This shock thickness can still be very thin to 
be resolved by the used grid and may lead to non-linear numerical instability. The 
standard techniques to obviate the difficulty are either to add explicit artificial dis- 
sipation or to employ an upwind biased scheme with implicit built-in dissipation. In 
the present effort the latter has been used via the MUSCL interpolation of variables. 

The conservative formulation for the Navier-Stokes equation in Cartesian 
coordinate system is given by, 

dq dF c d£c = 1 ,df v dG v . 

dt dx dy Re { dx dy j l “' j 

where, F c and G c are the convective terms in 2-D flow domain w r hile F v and G v are 
the dissipative or viscous terms in the flow domain. 

This Navier-Stokes equation has been solved for the flow variablesrdensity 
(p), velocities [u, u] r , pressure (p) and the total energy (E). 

The total energy(E) could be related to pressure as, 

V = (7 - 1)[£ - \ pW 2 + 1’ 2 )] 

The local speed of sound can be expressed as, 

P 


where 7 is the ratio of the specific heats C v and C v . 
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The unknown vector q is represented by, 

q = [ p , pit, pv, E] T 

To enable us to treat non-uniform grids or mapping a complex physical domain, we 
have worked in the transformed plane (£, 77) such that 

£ = £(s,y) 

77 = 77(3,2/) 

The main reason for working in the transformed plane is due to our strategy of using 
higher order spatial interpolation for accuracy. 

The governing equation (2.1) therefore transforms to, 

dq dF c 3G C _ 1 dF v . 

dt + d£ + d rj Re d£ + drj ’ (2 ' 2j 

where, 

Q = Q/J , 

F c = Uq + p { 0 €* U} T /J, 

G c = Vq + p[0 Tjx Tjy Vf/J, 

E v = {ixT+ yS)/J, 

G v = (VxT + fys)/ J, 

f=p(0 Tn r 2 i £7 i) t , 
s = p(0 r 12 r 22 ct 2 ) t 

with 

4 2 

Til = + TjxU-rj) g Vy^T])) 

r 2 1 = r 12 = f y it c + 7] y u v 4- f x i; ? + q x v v , 
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r 22 = g ~ g (6r«{ + *7 xWij), 

- UT n + V T12 + JTZijp + Vx{c 2 )ri\, 

AC 

a 2 = ur 2 i + vt 22 + JTTYjp ^( c2 )c + %( c2 )^] 
and the Jacobian of transformation and the contravariant component of the velocity 
in the transformed plane are given by. 


J ^xVy ^y^lx: 
Jj = + ^yV- 


V = T) X U + T)yV 


Here in the equation (2.2), Re is the Reynolds number, r)-s are the shear stress 
terms and a[s are the loss terms in the energy equation due to viscous terms. 
Enthalpy can be calculated as. 


H 


E + p 


Prandtl number is given by, 


Pr 




K 


where p is the viscosity which can be obtained by Sutherland’s [16] Law as, 

T Z/2 


p = 1.45 * 10 6 * 


T+ 110 


with T is the local temperature specified, C v is the specific heat coefficient at con- 
stant pressure and k is the coefficient of thermal conductivity. 

The Navier-Stokes equation can be reduced to Euler equation when Re — » 
oc, i.e. by retaining convective terms while neglecting the viscous terms. The Euler 
equation is, 

dq dF c dG c 

— -j -| = 

dt dl; drj 


(2.3) 
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The finite volume formulation for the equation (2.2) in terms of explicit time march- 
ing scheme for any arbitrary ii.j) th cell can be shown as, 

<7 n+1 - q n „ r - - l r . 

qJ. ~ J (Fcdr) - G c d£) + — j) ( F v drj - G v d£) = 0 '(2.4) 

which gives, 

^ ~ J 1 './ & dr] ~ (2.5) 

where the superscripts n and n + 1 show the present and future time step values of 
fluxes at any ( i,j) th cell. 

The time At represents the value of time step marching while S- is the 

v "' 

area of (z. j) th cell over which we are solving the Navier-Stokes equation. 

Here, the variables are referred to the cell center values. The integral rep- 
resents the flux integral performed over the cell interfaces. 

2.2.2 Nondimensionalization of the Quantities in Flow Equa- 
tions 

Nondimensionalization reduces the number of computations to a single parameter 
value. In arriving at equation (2.2) the following nondimensionalization has been 
used, 

Following scales have been used for nondimensionalization of the flow vari- 
ables: 

1. Reference Velocity: 

2. Reference free-stream properties: Poo , Uoo , Voc , £«,, M x , c an d Poc . 

3. Reference Specific heats and universal gas constant: C P0 , C VQ . Rq. 

4. Other reference parameters: L, Po ,Tq,k 0 
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For the computational cascade domain, nondimensionalised parameters can 
be calculated as given below. The asterix (*) quantities are the dimensional quanti- 
ties while nondimensionalised quantities are written without any super or subscript. 
The Tangential Velocity is calculated as. 


u = u m /U x 


The Normal Velocity is calculated as, 


The Density is given by, 


The Energy can be shown as. 


The Enthalpy is shown as, 


v = v*/U. 


P = P*/P= 


E = 


Pool'l 


H = 


H m 


Poet- i 


The Speed of Sound can be shown as, 


c = c*jc 0 


Using these nondimensionalized quantities we can express nondimensionalized vis- 
cosity by Sutherland's [16] law as, 


P = 


P* _ ,T* , 3 / 2 ^ Tp + 110 
Po _ To + 110 


U sing these parameters one can proceed for the calculation of Prandtl number and 
many other related parameters. 
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2.2.3 Nondimensionalizing the Inflow Boundary Conditions 

The computational domain with the specified boundaries has been shown in Figure 
(3.2). 

For the free-stream boundary condition the reference parameters are Poo.Uoo- Voo, 

and Eoq. 

For the inflow of the cascade domain, density is 
p = 1 as p* = poo. 

The velocities are, 
u=cos(0) and v=sin(0) 

as u* = IJ 0 o, v* = Uoo and 9 is the angle at which the flow enters the flow 

domain. 

The pressure can be shown as, 
p=l as p* = Poq. 

Hence, the free-stream energy can be calculated using nondimensionalized 
quantities as, 

E = + ip(u 2 + V 2 ) 

reduce to, 

1 1 

E b - 

7-12 

This shows that by nondimensionalization of the various parameters, we are making 
equation flexible for any large or small physical values. 
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2.3 Flux Vector Splitting Scheme 

2.3.1 Introduction 

All second-order central schemes generates oscillations which have to be damped by 
the addition of artificial dissipation terms to suppress the nonlinear instabilities. 

In smooth regions of the flow, where the flow variables can be considered as 
continuous, central schemes based on Taylor series expansion can be used with any 
order of accuracy. This is valid for supersonic flows also, where the apparent con- 
tradiction between the physical one-way propagation of waves and the symmetrical 
central difference schemes are direction independent, being resolved b}- considering 
the analytic continuation properties of smooth function. The validity of Taylor se- 
ries expansions express, indeed, a most remarkable property of continuous functions, 
namely that suffices the value of a continuous function in a single cell to be able to 
reconstruct the function in an increasingly large domain around that cell. In the 
limit, the complete local knowledge of its value and all its derivatives in a single cell, 
is equivalent to the knowledge of function everywhere. However, as soon as discon- 
tinuities appear this information is destroyed and more physical input is required in 
order to resolve the non-linear behaviour. 

The first level introduces only information on the sign of the eigenvalues, 
whereby the flux terms are split and discretized directionally according to the sign 
of the associated propagation speeds. This leads to the flux vector splitting methods, 
one of which we are using for our case of study in this thesis. 

2.3.2 Multi-Dimensional Flux Vector Splitting 

The flux vector splitting upwind methods can be extended to multi-dimensional 
flow by applying the one-dimensional splitting to each flux components separately 
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according to the sign of the associated eigenvalues. 

The Euler equation given in the section 2.2 for the transformed plane can 
be shown as. 


.4 = 


B = 



dq dF c dG c 
dt + d£ ' dq 

= 0 

(2.6) 

The convective terms from the equation (2.6) can be rewritten as, 



dF c dF c dq 
dE, dq dE, 

4 

(2.7) 

and 

dG c _ dG c dq _ 
dq dq dq 

4 

(2-8) 

where, 




> 

II 

4*' 

and B= 

dG c 

dq 


For the equations (2.7) and (2.8) the 

matrices A and B. which are the 

Jacobian Flux Matrices [3] are given as, 



^ 0 


& 

0 

- uU 

D + (2 

~ (7 - l)fzU 

(7 - 1 

^91 - vU 

c* 

& 

i 

cr 

i 

£ 

U ~ (2 - 7 )vf v 

(7 - 1 


7&- {uU + ^qj){j - 1) 

f (r -( v u+ftf)( 7-1) 

7 U 

and 




( 0 

Vx 

Vy 

0 

2 Y l VxQi - uV 

v + (2 - 7 ) 1177 * 

UVy ~ (7 - l)VxV 

(7- 

'-T-w! - vV 

vq x - (7 - 1 )q y u 

T ' 4- (2 - 7 )vq y 

(7- 

V b-l)V^-^£ 

?U. - (uV + 3f«?)(7 - 1) - (»V + \q \){ 7 - 1) 

jV 


\ 


The Jacobian matrices A and B have four eigenvalues which can be repre- 
sented in terms of right and left eigenvectors (R and S) as, 
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\ A = RAR~ l = 


and 


0 

0 


0 

U 

0 

0 


0 

0 

u 

0 


0 

0 

0 


u + cy/g+% 


(2.9) 


A b = SBS~ l 


0 

0 

0 


0 

V 

0 

0 


0 

0 

V 


0 

0 

0 


0 V + Cy/fjf+rj* 


( 2 . 10 ) 


where. R and i? _1 are the right and left eigenvectors whereas S and 5 -1 are 
the corresponding eigenvectors for B-matrix. 

These eigenvectors can be shown as, 


( 


R = 


1 

V. 

V - iyC 


0 


-tv 


1 

u 


i \ 

XL A C 
V-iyC 


\ H - © C c -iyU + i x v \{u 2 + V 2 ) H + Q^c ) 


( 2 . 11 ) 


R~ l 


1 

2 


/ (6 2 + ^) + ~(b 1 v + ^) h ^ 

2(£yU-£xV) -2 4 24 0 

2(1 — fe 2 ) 2biU 2 b } v -2bi 

-(M-&) 6, j 


( 2 . 12 ) 


and similarly, 
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S~ l 
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% 
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rjxC 
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VyC 
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~Vx 

r -f 

TjyC 



U- 

- Q v c 

W 2 + V 2 ) 

TlyU - 

rjxv H d 

QyC J 


/ 

(&2 + 

C ' 

—{biU + Vf) -(biv—^-) 

b 1 

1 


2 ( 1 - 

62 ) 

2b\U 

2b\V 

—26], 

2 


2 (jf*u 

- Jfc«) 

2 7 ?y 

—2 

Vx 

0 


U 6 ;- 

© 2 .\ 

C J 

~(M - 

*) 

61 


\ 


J 


where, 


(2-13) 


(2.14) 


h = + v 2 ) * b u 

c hi 

^ = /cTZT 5 : 

+ f y v 


and 


- _ Vx 

V^I + ’ 
% 

Vy = 7=f= p 

V 7 ?! + ^y 

©r, = Tfi'U + TfyU 
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Hence, an upwind formulation can be obtained with the Jacobian matrices as, 



A ± = RA^R - 1 

(2.15) 

and 

B ± = SA±S-' 

(2.16) 

with the Jacobian matrix property of, 



A = A + + A~ 

(2.17) 

and 

i 

1 

+ 

(2.18) 

similarly, 

B = B + + B~ 

(2.13) 

and 

\B\=B + -B~ 

(2.20) 


The fluxes associated with these split Jacobian matrices are obtained from the re- 
markable property of homogeneity of the flux vector F c (q). The function F c (q ) is a 
homogeneous function of degree one of q. 

Hence the flux splitting can be defined as, 

F c = Aq = F* + F~ 


and 


G c = Bq = Gt + G- 
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where 


f+ = .vr, 

F- = A-q- 

(2.21) 

Gt = B + q+, 

G- c = B-q- 

(2.22) 


The splitting of flux vectors as shown in equations (2.21) and (2.22) are directed by 
the property of information propagation that can be achieved by upwind schemes. 

The equations (2.21) and (2.22) show. A + as the Jacobian matrix corre- 
sponding to the disturbances propagating from left to right and A~ corresponds 
to the disturbances propagating from right to left. Similarly. B + corresponds to 
disturbances propagating from bottom to top while B~ from top to bottom. 

In this form of flux, vector splitting q * is the flow property just to the left 
of the surface considered while q~ is just to the right of the surface considered for 
equation (2.21). while for the equation (2.22). q * is just to the top of the surface 
considered and q~ is to the bottom of the surface considered. 

In actual computer coding for our convenience, we are assuming these split 
flux quantities as. q L . q R for left and right side while qr- q R for top and bottom side 
of the surface considered. 

The Figure (2.1) shows the flux vector splitting for ( i,j) th cell for the left 
and bottom surface of the cell considered. 

The eigenvalue matrices shown by equation (2.9) and (2.10) decides the flow 
nature i.e. subsonic or supersonic at an}' cell in physical domain. The flux splitting 
technique by Steger and Warming [4] that has been used here is explained in next 
subsection. 
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Figure 2.1: Flux vector Splitting for {i.j) th cell 


2.3.3 Steger and Warming Flux Vector Splitting 

From the equation (2.21) and (2.22), we can calculate the eigenvalues (A) for 
A~ and B~ . B~ matrices by the method suggested by Steger and Warming [4] as. 


A + _ A++ 1 A ‘ 1 

2 

(2.23) 

A"- 1 A" 1 

(2.24) 

A '= 2 


For each of the eigenvalues of A 
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where, 


Af 


Af 

0 

0 

0 


0 

Aj 

0 

0 


0 

0 


a 3 


0 

0 

0 

Ai 


(2.25) 


and 


Af 


(£ +Cy/(* + £g)± [ {U 4- cyg + Q) | 
2 


Af = 


Af = 


U± 1 U I 
2 

U± \ U I 


Af 


_ ~ C \/^i + | (U - Cy/g + £)) 1 


where, 




Af 

0 

0 

0 


Af 

0 

0 


0 

0 


a 3 


0 

0 

0 

Af 


Af = 


(V + cy/rg + T) 2 y )± | (V + Cy/irg + r% 


(2.26) 


Af = ^ 


±_Z±Jil 


a 3 = 


Af 


(F - Cy/rg'+ig)± | (V - c^/rg + rg) | 


2 
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Here. A : has only positive eigenvalues. A only negative values and such 

that'. 


A = A + + A" , |A|=A + -A“ 


2.3.4 Properties of Split Flux Vectors 

The most important concept of flux splitting is that, it is totally dependent on the 
fact that the fluxes are homogeneous function of degree one in q. One can not 
therefore directly apply this approach to a general scalar flux function f(q). Since 
the only homogeneous scalar flux function of degree one is the linear case / = aq. 

As shown in equations (2.21) and (2.22) splitting of Jacobian matrices give. 


A = .4 + + .4“ 


(A?7 


7) 


and 

B = B + + B~ (2.28) 

with 

F c = Aq — F c + + F c ~ = A + q A~q (2.29) 

G c — Bq — G c + G c = B ' q ~r B q (2.30) 


and 


d _I± 

dq 

dG c 


A = 


dF+ dF, 


dq 


+ 


*r = A+ + A- 
dq 


= B = =B + tB' 


dq dq dq 

but one does not generally have the equality of the split Jacobian i.e. 


1.31) 


(2.32) 


djt 

dq 





(2.33) 
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or 


aG ‘ , B+j 


dG: 


^B- 


(2.34) 


8q ' ~ 7 dq 

These matrices do not have the same set of eigenvalues. The eigenvalues of 
„4 + and A~ as well as B + and B~ arc the positive and negative values of A and 13, 
that is A* and , but this is not true for and 0G * 


dq 


The non-equality of and with A ± and B ± has very important 
relation with regard to the definition of upwind schemes based on the concept of 
flux vector splitting. Therefore, the conservative equation in split form can be shown 


as, 


dq , dF t + , OF- , d&t , dG: 


+ 


+ 


+ 


+ 


= 0 


dt d£ ’ df; ' dr] drj 
which in terms of quasi-linear form can be expressed as, 


(2.35) 


dq dF+dq dF~ dq dG+ dq dG: dq n , n 

i + wi + -el-i + -ati + ^ti =0 (2 - 36) 

Alternatively, if we write the quasi-linear form first and then split the Ja- 
cobians then, 


df I , , B d V - d( l , A + d V , A - d Q , B +d<! , B -dq _ () 

Tt +A TG at +i4 + di d~G dj~ 


(2.37) 


The two formulations are not identical, as a consequence of equations (2.34) 
and (2.33) showing that the two operations of splitting fluxes and introducing the 
quasi-linearization do not commute. Equation (2.37) keeps the same physical char- 
acteristics which are separated accordingly to their sign while equation (2.36) does 
not represent a decomposition of the physical characteristics, but is the only lin- 
earization of consistent with the conservative form of equation (2.3G). 
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2.4 Interpolation Scheme 

2.4.1 Introduction 

The motivation behind the interpolation schemes is the hope that the introduction 
of physical propagation properties in the discretization will prevent the generation 
of oscillations in the numerical solutions. This is only partly fulfilled in the sense 
that for nonlinear equations, such as Euler equation oscillation free results can be 
obtained for weak stationary discontinuities. However, this is not a general property, 
since it can be shown theoretically that the linear second-order upwind schemes al- 
ways generate oscillations. The concept of monotonicity was introduced by Godunov 
(1959) to achieve the oscillation-free solution. 

As a piecewise numerical approximation is equivalent to a first-order spa- 
tial discretization, a linear approximation on each cell gives second-order spatial 
discretization: while a quadratic representation on each cell leads to third-order 
spatial discretization. Here in this thesis, we have used third order MUSCL type 
of spatial discretization which has built-in dissipation for solving full Navier-Stokes 
equation. 

2.4.2 High-Order Upwind Scheme- MUSCL approach 

Second-order spatial accuracy can be achieved by introducing more upwind points in 
the schemes. The state variables at the cell interface or at the cell centre are obtained 
from an extrapolation between the neighbouring cell averages. This method for the 
generation of second-order upwind schemes via variable extrapolation is referred as 
the MUSCL approach. MUSCL acronym stands for, Monotone Upstream- centred 
Schemes for Conservation Laws, as developed by Van Leer (1979). The time inte- 
gration methods based on a separate time and space discretisation , such as multi- 
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Figure 2.2: Flux Estimation for the cell 

stage Runge-Kutta techniques or linear multistage methods (implicit schemes), can 
be directly applied to the MUSCL scheme. 

The representative cell for the flux estimation is shown in Figure (2.2). 
where, q is the flux quantity and the superscripts L and R refers to the left 
and right sides of a cell surface while, T and B are the top and bottom sides of the 
cell surface considered. 

Higher-order upwinding scheme with the generalised formulas can be written 
as, 

for the left hand side of the cell, 

q L i-i! 2 j = Qi-ij + |[(1 - *)(&- ij - Qi-ij) + (1 + «)(<7t,j - Qi- ij)] (2-38) 

Q R i-i/2j ~ QiJ - |[(1 + _ Qi-u) + (! - - Qij)] (2-39) 

Similarly for the right hand side of the ( i,j) th cell we have, 

Q L i-i/2j = Qi j + | [(1 _ K )(QiJ ~ Qi- ij) + (1 + K )(Qi+ij ~ Qij)] (2.40) 
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Figure 2.3: Third order interpolation within the cells 

1 / 2 .J = Qi+ij ~ |[(1 + k)(?.-+1j - Qij ) + (1 - «)(ft+2j “ ft+ij)] (241) 
While, for top and bottom side, generalised interpolation formula is given by. 

Q T ij- 1/2 — Qij - 1 + ^[(1 — K ){Qij-l ~ Qij- 2 ) + (1 d- K ){Qij ~ Qij- 1)] (242) 

q B ij- 1/2 = 9.-J - |[(1 + - Qij- 1 ) + (1 - - 9ij)] (243) 

Q^ij- 1/2 = T ~[(1 — K )(Qij ~ Qij-i) T (lt K )(Qij+i 9t.j)] (2.44) 

Q B ij+i : 2 = Qij-i - 1 — ^[(1 + - 9ij) + (1 - K)(qij + 2 - Qij-ri)} (245) 

where, the introduction of parameter e, such that e=0 is for a first order upwinding 
scheme and e=l is for a higher order upwinding schemes. While k defines the nature 
of the scheme given by Hirsch [9] as, 

k — — 1 => One sided linear interpolation scheme 
k — 0 = 4 - Linear interpolation scheme 
k = 1 =4> Central differencing scheme 
k = 1/2 =£- QUICK scheme 
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and finally, 

k = 1/3 => MUSCL interpolation scheme 

In this thesis, we have used e=l and k— 1/3 for third order MUSCL 


inter- 


polation scheme. 

From equation (2.38) to (2.45) substituting e=l and k—1/Z 


we get, 


„L , 9:j 9i-lj , 9i-lj 9i-2J 

9 i— 1/2 j — 9:-lj 1 o "h R 


Qi—lj QiJ , 9ij 9i+lj 

9 i— 1/2 j - 9i,j + 3 h g 

„L Qi+lj ~ Qi,j , Qi,j ~ Qi-lJ 

9 i+l/2 j - Qij + 3 i g 


„K . 9:j Qi+l,j , 9i+lJ 9i+2j 

9 X+1/2J — 9t+i,j V 3 1 g 

; 9tj ~ Qi,j—1 Qi,j- 1 ~ 9ij-2 
f-i i x 1 r 


: , 

9 jj— 1/2 — 9ij'-l V ~ T 

„B 


3 ' 6 


„B , 9ij-l 9i,i , 9i,i Qi,j+1 

9 ij— 1/2 - 9ij + 3 r - 

T _ _ , 9i.j'+l ~ 9iJ , 9iJ ~ Qij-1 

9 tj-fi/2 Qij * 2 * g 

-B _ „ , 9ij ~ 9i ,j-rl 9i,j+l ~ 9i,j+2 

9 ij-fl/2 9x J+l V - o 


(2.46) 

(2.47) 

(2.48) 

(2.49) 

(2.50) 

(2.51) 

(2.52) 


2.4.3 Free Stream Preservation for the Flux Calculations 

The Euler equation for 2-D compressible, unsteady and uniform flow is rewritten 
for the transform plane as, 


dq dF c dG c 
dt d£ dq 

For an uniform flow w r e have, 


(2.54) 


q[ — k \ , ?2 — &2) 93 — ^3j 94 — 


(2.55) 
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A [hoi 0: Poc^oc: Poc^'oc) -E"oc] 

where, oc is the free stream value of the flow variable. 
From equations (2.55) and (2.56) we have. 


(2.56) 


Qi — hi/ J.qo — A'2 / J. <73 = A 3 / J, <74 — /c 4 / J 

where for a single {i,j) th cell Jacobian(J) will remain constant for all four 
quantities. 

The convective terms from equation (2.54) can be rearranged as. 


similarlv. 


dGc 

dr] 


= 4 ®i= 4-(-) 
dd, dCJ j 

(2.57) 

= = 

dr] dr] J 

(2.58) 

dq d , k 

dt ~ dvr 

(2.59) 


From equations (2.57), (2.58) and (2.59) we can rearrange the 2-D Euler equation 
as, 

r) Jr r) Jc r) h 

(2.60) 


d k d .k <9 ,A. 

S ( j ) + -Vi )+ % ( j ) = 0 


But 57(7) =0 as J ^ J(t) 

Thus the equation (2.60) becomes, 


d ,k. _ d .k 

A dZ^ + B dvty ~ 


(2.61) 


In case of 1-D flow' as shown in Figure (2.4), 


, _ (lAQi+l ~ = n 

d£'J ll 2A<£ 


(2.62) 
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i-2 i-1 i i+1 i+2 

Figure 2.4: Cell distribution in 1-D Flow 


which means. 

Ji+i = Ji—i ( 2 . 63 ) 

The equation (2.62) tells us that, for uniform flow the jacobians are in checkerboard 
mode i.e. the alternative cell Jacobians must be same for the free stream preservation 
which is difficult to have in practice. 

Similarly for 2-D flow equation we have checkerboard mode in as well 
as in rj direction. This shows that the grid plays major role in compressible flow 
calculations. 

2.5 Time Integration Scheme 

2.5.1 Introduction 

To obtain a steady state solution for the flows at high Reynolds number, the time- 
step restriction becomes a severe impediment to overall efficiency. This time step 
restriction appears as the Courants(CFL) number criteria for the flow calculations. 
This difficulty can be overcome by introducing Runge-Kutta [1] time marching 
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schemes which allow larger CFL numbers. For unsteady flow problems, it would 
be appropriate to use higher order Runge-Kutta schemes which are stable for CFL 
number C < 2\/2 if the dissipative terms are small enough. 

The Rational Runge-Kutta scheme(RRK) for Euler as well as Navier-stokes 
solver has been developed by Morinishi and Satofuka [10] . For the flow domain, the 
accuracy of the solution can be improved by increasing the number of grid points. 
Increasing the number, on the other hand, may reduce the efficiency of the solution 
considerably. Since a smaller time step requires a greater number of iterative steps in 
addition to increase in operations counts per step. This difficulty has been overcome 
by Morinishi and Satofuka [10] by introducing RRK time stepping. 

In this thesis work, the two stage RRK scheme has been used as a time- 
stepping scheme. With this scheme the solution is fully explicit as well as robust. 

2.5.2 RRK Formulation 

The Xavier-Stokes equation discretized by Finite Volume principle, as shown in 
section 2.2. can be rewritten in simplest possible way as, 

§ = <3(«) (2-64) 

For the integration of this ordinary differential equation, the two stage RRK 
scheme [10] is used. The two-stage RRK scheme can be written as, 

9i = A iQ(q n ) (2-65) 

92 = A tQ(q n + c 2 gi ) 

93 = b\9i + bog2 

[2gi(giWb) ~ 93(.9i-9i)} 

^ ? {93-93) 


( 2 . 66 ) 

(2.67) 

( 2 . 68 ) 
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where, the superscript n denotes the index of time steps and (5,. g 3 ) denotes 
the scalar product of two vectors g { and §j. The coefficient 6 1; b 2 and c 2 must satisfy 
the relations: 

b\ — b 2 = 1 and & 2 c 2 < —0.5 

For an efficient second order scheme. 

6i=2, b 2 ~- 1 and c 2 =0.5 which have been used here. 

The time step At is determined locally so that the local CFL number be- 
comes constant for each grid point. The critical time step for each coordinate direc- 
tion is obtained as follows: 


(At^)ij — [- 


1.0 


(At n )ij — [' 


u | +C^ T 2 + £, 2 

1.0 


1 . . 
J 1 J 


(2.69) 


(2.70) 


I V 1 +t\4 2 + %- 

The critical time steps are used in the estimation of dissipative fluxes. The net time 
step for each grid point is determined from these critical time steps as, 


Atij = C 0 min(At A t T] ) i i 


(2.71) 


Here C a is the local Courant number and c is local speed of sound. 

The viscous terms of the Navier-Stokes equation and the dissipative terms 
are calculated only in the first stage of RRK scheme and frozen in the second stage. 
The boundary conditions are called after even- stage. 

Finally, we could write the equation (2.68) in the form of, 


(7 n+1 = < 7 ” ■ 4- f • ■ 
9i j ~ 9jj ~ 




(2.72) 


where, 

- [2gi (gi , 93 ) - £3 (9i ; 9i j] 
r,J ' (as, gz) 


(2.73) 
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This fj j term called as a residual term in two-stage RRK scheme. For the fast 
acceleration of the convergence, residual terms play major role, which could be 
handled by various accelerating techniques, which are not used in this thesis. 

2.6 Boundary Conditions 

Proper boundary conditions enforcement is the most important part in the compu- 
tation. Particularly, for internal transonic flows, the boundary conditions are critical 
for the correct solution to be computed. 

Here, in this section, for a given cascade domain, four boundary conditions 
have been discussed. As shown in section 3.3.2. Dirichlet boundary condition has 
been enforced for density and momentum quantities at inflow boundary. For our 
case of physical domain, the following are the broad categories of the boundary 
conditions as shown in Figure (2.5). 

(i) Periodic boundary conditions at AB, EF and similarly on CD and GH. 

(ii) Solid wall boundary conditions at BC and FG. 

(iii) Far-field boundary conditions 

a. Inflow boundary conditions on AE. 

b. Outflow boundary' conditions on DH. 

2.6.1 Periodic boundary condition 

As shown in Figure (2.5) , the boundaries AB, EF and CD, GH are the boundaries 
where periodic boundary' condition is implemented. Since, these boundaries are in 
common free flow boundaries to the other rows of turbine blades, the flux exchange 
at the cell j—1 will be same as at dummy cell j=m, which will be the first cell in 77 
direction for second row of the turbine blades. 
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Figure 2.5: Periodic Boundary Zone 



2.6.2 Solid wall boundary condition 

The Figure (2.5) shows, boundaries BC and FG as the solid wall boundaries in trans- 
formed plane. Solid boundaries are the actual blade zones. This zone is being formed 
with NACA0012 and DCA aerofoils for two different cases of study undertaken here. 

For calculations of C p values at upper and lower surface of aerofoils in Euler 
as well as in Navier-Stokes code, the boundary condition at solid wall is important. 

On a solid wall, while computing inviscid flow in Euler code, it is necce- 
sory to satisfy the zero normal velocity boundary condition as shown by Ravichan- 
dran [11]. Similarly, for a viscous flow calculations both the tangential as well as 
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normal velocities should be forced to zero. 

The contravariant component of tangential velocity near the '.vail is given 
as. 


U = + v£ y (2.74) 

while normal velocity is given as. 

V' = ur] x + vrjy (2.75) 

Euler solution boundary conditions: 

Since, present scheme is of cell center type, the solid wall forms a cell interface 
between j — l and j = 0 latter being the first of the two dummy cell rows introduced. 
A simple way to satisfy the zero normal velocity condition in inviscid flow on the 
cell interface is to use the reflection principle. This principle says that the normal 
velocity should be antisymmetric i.e. 

lj=o = — 1 j=l 

while other quantities (p, U, E) will remain symmetric for all dummy cells along the 
solid wall interface. The Figure (2.6) shows the reflection principle graphically. 
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Figure 2.6: Solid Boundary condition for Euler calculations 
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With the above mentioned solid wall boundary condition, we get the values 
for Uj-o and Vj - 0 as. 


„ Lj=lVy io • iji=lCy |o 

U j = 0 ~ r 

Jo 

Lj- \T) X |o -fl'j = x^ r |o 

V 3 = o 7 

Jo 

where. J is jacobian matrix and relation between rj x . £ y and rj y is. 

(si)j= o = (fi)j=l 

(Cy)j=0 = (£y)j=l 
(Pz)j=0 = {Vx)j = 1 
{ r jx)j = 0 = (' T Jx)j=l 

and 

Jj=0 = A/=l 


(2.76) 


(2.77) 


Navier-Stokes solution boundary conditions: 

Similarly, for Navier-Stokes solution both the contravariant velocity components 
should be antisymmetric across the solid surface for validation of reflection principle. 
As the flow is viscous in nature near the solid surface we treat tangential velocity 
component also to be antisymmetric, i.e. 


Uj = o — —Uj = i 


and 

Vj=o = —Vj = i 

While other quantities (p and E) are still retained symmetric across the solid surface 
boundary as shown in Figure (2.7). With the above mentioned boundary conditions 
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-V 


Figure 2.7: Solid Boundary condition for Xavier-Stokes calculations 


the equation (2.76) and (2.77) get modified for Xavier-Stokes boundary conditions 
as. 


L’j^Tjy |o +^j=l£y |o 

U J=0 - J 

o'O 

(2.78) 

&j~iT] x |o j=l £rr |o 

Jo 

(2.79) 


with all other constraints remaining same as of the Euler solver. 

Requirement of the dummy cells is dependent on the order of the upwinding 
scheme used near the cells interface. For a third order MUSCL interpolation, we 
required two dummy cells while for a first order upwinding just one. 

Temperature Variation across the Solid Surface: 

In case of Xavier-Stokes calculations temperature at solid boundaries is also the 
important boundary condition needed. 

For our case of study, we are forcing an adiabatic condition on the surface. 
According to the adiabatic condition, the temperature gradient along the normal 
direction to the solid surface will be zero. As shown in Figure (2.8) the normal 
vector to the solid surface can calculated as, 


n = isin{9} — jcos{6} 


(2.80) 
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Figure 2.8: Normal Vector profile for the Solid Boundary Surface 


where, 


sin{6} = Ay/ As,cos{6} = Ax/ As 


The temparature gradient for adiabatic condition can be calculated as, 


dT 

dn 


VT ■ ft = 0 


since |^ = 0 is for adiabatic condition, which gives, 


I - ms{e} % = ° 


where, 

dT _ dT dT 
dx ~ dC x+ d V Vx 
dT _ dT _ dT 

Jy~ dC y 1 dy ny 

Since we know the Tf at any point in solid boundary zone we can find out 
equations (2.82) to (2.84) as, 

_ T^sinjO} - g y cog{g» 

17 {Vx sin{9} - rj y cos{6}) 


(2.81) 

(2.82) 

(2.83) 

(2.84) 
T v from 


(2.85) 
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2.6.3 Far-field boundary conditions 

Far-field boundary conditions for the present case are inflow and outflow boundary 
conditions. For inflow as well as at outflow boundary conditions. ID Riemann 
Invariant [7j conditions are used. 

Inflow Riemann Invariants 

For the subsonic inlet condition, three Riemann Invariants are to be evaluated ac- 
cording to the free-stream value at upstream as they carry the information from 
outside into the flow domain while the fourth quantity is extrapolated from inside 
to the inflow boundary. In case of the transonic inlet condition the three flow vari- 
ables (p. pu. pv) are evaluated from free-stream values and the fourth variable Energy 
(E) is calculated from three free-stream variables. 

Inflow Riemann Invariants are calculated as [7], 


Ri = (W (2.86) 

P 7 

R 2 = Qtoc ( 2 - 87 ) 

Rz = (?„ + ~^r)o o (2.88) 

7-1 

i?4 = (?n - ~~~r)e (2.89) 

7-1 


where subscript oo refers to as free-stream value, e is internal flow domain value, t 
is the tangential flow to the inflow boundary' and n is the normal direction of the 
flow to the inflow boundary. 

With the equation (2.88) and equation (2.89) we get the local normal ve- 
locity and speed of sound at inflow boundary as, 

Qn = ^(-^3 + Ri) 


(2.90) 
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c= 2_i( jR 3_ i?4 ) (2.91) 

4 

With these quantities we can obtain Energy at the inflow boundary. 

Outflow Riemann Invariants 

On a similar basis, we can calculate the outflow Riemann Invariants, where we have 
pressure as a fixed quantity. With fixed pressure at outlet, we can easily find out the 
Energy at outflow while other three quantities have to be extrapolated from inside 
the domain. 

Outflow Riemann Invariants could be shown as. 


P — Pspecified 

(2.92) 

R ■ - (y>« 

(2.93) 

Rt. = Qt e 

(2.94) 

_ , 2c . 

#3 = (<?n + . )e 

7-1 

(2.95) 


These Riemann Invariants give the correct interpretation of the flow at outflow 
boundary. 

2.7 Spectral Analysis of FVS Scheme 

2.7.1 Introduction 

It is well known for incompressible flows [14], any discretization scheme other than 
the one using spectral method introduces phase error and numerical dissipation. For 
high Reynolds number even spectral methods require numerical dissipation. While 
numerical diss ipation is necessary for computing flows at high Reynolds number, it is 
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the phase error which determines the eflieacy of the numerical method. The source 
of such errors are due to the treatment of the nonlinear convection terms in the 
conservation equations. We are not aware of similar analysis for compressible flows. 
In this section we extend the ideas in Sengupta & Sengupta [14], for the spectral 
analysis of the flux vector splitting schemes that are used in finite volume methods 
for flic solution of compressible flows. Although the analysis is for the convection 
terms alone, the results would be valid for full Navier-Stokes solution. The analysis 
assumes greater significance since we use a higher order interpolation for the split 
fluxes. 


2.7.2 Spectral Analysis 

We analyse the schemes with respect to equation(2.2) given before. Here the inviscid 
fluxes are written in terms of the flux Jacobians A and B. 


dF c dq 

A d$ 

(2.96) 

dG c _ B d Q 

(2.97) 

dq dq 

The unknowns q are written using Fourier Spectral representation 

as given by, 

q(thV) = [ qi{k,q)e iKtil dk 

(2.98) 


where & = (l — 1)A£ 

Here we will only show the analysis with respect to one variable £, similar 
ideas can be used for the other direction as well. 

After applying the finite volume technique, we actually solve the euqa- 
tion(2.4), where one of the fluxes are written as, 


F c = Aq 


(2.99) 
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B 


<l,m) 


D 


Figure 2.9: (l, m) th cell 

For a given (l, m) th cell as shown below, let us illustrate the procedure 
above flux given by equation(2.99) is represented by, 

Aq = A + q + + A~q~ 


for the side AC. 

In equation (2. 100) the left hand side is evaluated as, 

Aq=Aj ft c W*-3/2)A t dk 

On the right hand side, the flux Jacobians are given by, 

A = A + +A~, B = B + +B~ 

with 

A + = ra a + r~\ A~ = rAa-r- 1 

The various quantities are already defined in equation(2.11) to (2.16). 
Similarly the unknown q are split according to 

q = q + + (f 

with the split unknowns given as in equation (-2.42) and (2.43). 


when the 

( 2 . 100 ) 

( 2 . 101 ) 

( 2 . 102 ) 

(2.103) 

(2.104) 
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Therefore the equation(2.100) can be rewritten as, 

+(fiA„-J{- 1 )[5, im - { j, +lm _ km}] 

(2.105) 

where we have used the general representation for q + and q~ as given in Hirsch [9] 
Since A = RA a R~ 1 , the above equation can be also written as, 


ActualPhase 
EquivalentFV Phase 


A r ~»*A{ e(l — k) . -»«A£ -3i>cA; e(l + k) , i*At -i*A{ 

9 I J <) /i O ll ' ^ 1 « o __ g 2 11 


A /i + At 


4 [e"T“ + ^-A{ e T“ - e-^} + AAJ_A{<r 


A, r e(l + k) r xk efl — /c) r 3i*Ag ikA£... 

+ A - + X+ I e 1 1 “ e ’ > ^ ’-e'}} (2-106) 


Ideally the left hand side should be equal to unity. The used splitting shows 
that the right hand side is a combination of a real and an imaginary quantity. The 
departure of the real part from unity represents the phase error or wave number dis- 
tortion. Physically such distortion gives rise to numerical dispersion. The imaginary 
part contributes to the added numerical dissipation, which is absolutely necessary 
for numerical stability. 

In equation(2.106), various schemes can be simultenuosly analysed depend- 
ing on the value of e and k, which are discussed in details in section (2.5.2). For 
example, for the MUSCL interpolation, e = 1 and k = 1/3 as shown in equation 


(2.50). 

A+ 

Furthermore in equation(2.106) the factor is always positive if A is 

positive and greater than one. If this factor is denoted by a(> 1.0), the other factor 
in equation(2.106) will be given by (1 — a). 
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In Figures (2.10) and (2.11). the real and imaginary part of equation (2.106) 
have been shown for all a values i.e. irrespective of the value of the a, the real and 
imaginary solutions for different finite volume schemes remain identical. 

Certain features emerge from this figure. Firstly the abscissa is plotted for 
(kA£) up to a maximum value determined by the Nvquist criteria. Second feature 
is that the first order upwinding scheme(e = 0) and the central differencing(K = 
l.e = 1) show identical phase error, but the distortion is rather high. Notice. that 
the first order upwinding is nothing but the second order central differencing scheme 
with implicit addition of numerical dissipation-as can be seen from the imaginary 
Figure (2.11). Thirdly some of the other schemes including the MUSCL interpolation 
scheme shows some overshoots in the plotted range. While this is an undesirable 
attribute of any scheme, for the used MUSCL scheme(e = 1,k = 1/3) , this is 
minimum. The maximum overshoot is given in Table (2.1). 


i 1 

K j € 

' 

Scheme 

Max. Overshoot 

0 

i 

Linear Interpolation 

10% 

1 1 

Central Differencing 

no overshoot 

1/2 1 1 

QUICK 

no overshoot 

1/3 | 1 

MUSCL interpolation 

1% 

1 

- 1 

i 

One sided linear interpolation 

40% 

i Any 

0 

First order upwinding 

no overshoot 


Table 2.1: Phase error for different Interpolation schemes 


It is evident that the one-sided linear interpolation (e = 1 ,k = -1), is 
practically useless as it has tremendously high built-in numerical dissipation. At 
the same time the n um erical dissipation of this scheme is also the highest among the 



2.7 Spectral Analysis of FVS Scheme 


45 


plotted schemes at the high wave number range. Fourthly, when one compares the 
MUSCL interpolation scheme with the first order tipwinding scheme, one notices 
that the latter damps across all wave numbers, while the MUSCL scheme only 
damps the higher wave numbers.. High Reynolds number calculations requires the 
large eddies to be simulated without dispersion, while the numerical instability be 
removed by damping the high wave number components. Also notice that the linear 
interpolation (re = 1, e — 1), has phase error overshoot intermediate between MUSCL 
and one-sided linear interpolation. This case also has numerical dissipation that 
is in between the first order upwind scheme and MUSCL interpolation scheme. 
Finally, the QUICK scheme(re = 1/2, t = 1), seems to be optimal i.e. this scheme 
does not have any overshoot, at the same time has lower dissipation at lower wave 
number range. However if the added dissipation at the higher wave number range 
is adequate or not, to suppress nonlinear instabilities, can only be investigated with 
respect to a realistic simulation at high Reynolds numbers calibrated with respect 
to experimental results. 
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Figure 2.10: Phase error comparison for different schemes 



Figure 2.11: Numerical dissipation comparison for different schemes 





Chapter 3 


Grid Generation 


3.1 Introduction 

Grid generation is of importance to solve any flow problem for a given physical 
system involving complex geometries. Grids that are used can be classified into two 
major categories, 

(i) Structured grid, (ii) Unstructured grid. 

Structured grids are basically used for solving flow equations based on Finite 
Difference or Finite Volume principles, while unstructured grids are best suited 
for Finite Element applications where there is no definite structure is maintained 
between nodal points. In this work we are concerned with the former. 

Grid generation is a technique to establish the correspondence between point 
(x,y) in the Cartesian coordinate system and point (£, 77 ) in the transformed compu- 
tational domain. The problem of grid generation can be posed as a boundary value 
problem. The Figure (3.1) shows the boundary value problem has been based with 
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Dirichlet boundary conditions. In defining relationship between points in physical 
and computational domain i.e. x=x(f,77), y=y($.r/), it is necessary that there be a 
one-to-one correspondence. It would be unacceptable for a single point in physical 
domain to map into two or more points in computational domain and vice-versa. 
After establishing the mapping i.e. x=x(f, 77), y=y(£. 77). the requirement of one-to- 
one mapping can be determined by evaluating the determinant of transformation 
matrix. J . also known as Jacobian. For one-to-one mapping | J | must be finite and 
nonzero at each and every point in the grid domain. 

In this work, the flow equation based on the Finite Volume principle is 
solved using structured mesh with domain of (193x179) grid points for NACA0012 
as well as DCA cascade domain as shown in Figure (3.3) to (3.5). The grid for the 
cascade problem has been generated by an algebraic grid generation technique and 
by elliptic grid generator. 

The grid generated for NACA0012 aerofoil with (193x179) grid points by 
an algebraic technique is shown in Figure (3.3) in multiple domains. However such 
multiple domain grids will have discontinuous metrics. In this respect an elliptic 
grid generation method gives good grid smoothness and can be easily controlled for 
efficient grid spacing. 

The grid generated for the cascade by an elliptic grid generation technique 
is shown in Figure (3.4). An elliptic grid has been found more suitable for solving 
Euler and Navier-Stokes equations. 

3.2 Algebraic Grid generation 

3.2.1 Introduction 

The structured grid generation can be divided in to following categories, 
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(i) Algebraic method 

(ii) Conformal mapping 

(iii) Method using PDE’s 

An algebraic method is efficient and can be applied to any 2D or 3D physical 
systems. This kind of mapping extrapolates the boundary data to generate interior 
grids. The major requirement is that the generated grids should be well conditioned 
i.e. smoothly varying, close to orthogonal. By means of some explicit grid con- 
trol technique, we can have a desired grid spacing at critical points or zones in the 
domain. However, the slope discontinuities present at the boundaries tend to prop- 
agate to the interior. In such cases, the grid suffers from lack of smoothness. The 
distribution of the points along the boundary of the domain is handled effectively 
by defining normalized one-dimensional stretching functions along the boundary 
segments. One dimensional stretching function is described in next section. 

3.2.2 One-Dimensional Stretching Functions 

In an algebraic grid generation method, the transformation could be done analyti- 
cally where the boundaries are regular. However, an interpolation technique has to 
be used for an irregular complex boundaries. 

For the cascade geometry, we have treated whole domain as three blocks 
structure as shown in Figure (3.2). 

(a) An upstream zone(I), 

(b) Solid boundary zone (II), 

(c) A downstream zone (III). 

We have generated the grid algebraically in these multiple blocks. As the 
algebraic grid generation process is very fast, the mutiblock nature can be adapted 
for its suitability in solving the governing equation by domain decomposition, which 
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could be used in future for parallel computing. 

Here, the grids are generated in domains using normalised one dimensional 
stretching functions, ratio analysis and bidimensional explicit interpolations. For 
the cascade, we consider the whole domain as union of three non-overlapping do- 
mains and grids are generated in each domain separately. Let us briefly discuss the 
stretching function that is used for distributing non-uniform points on the domain 
boundaries for the first block. 

A normalised independent variable is defined as, 

* = (3-1) 

so that 0< <1 as £ a < £ < ft; 

where a and b are the first and last points along the boundaries of the 

domain. 


A stretching function is defined by Eiseman [1] as, 


a = Fr + (l-P)[l- 


tanh[Q( 1 — £*)] 


(3.2) 


tanhQ 

where P and Q as grid point control parameters. P effectively provides the slope 
of the distribution close to £*=0 ; Q is called a damping factor and controls the 
variation from linear dependence. Once s is obtained, it is used to specify the 
distribution of x and y. For example, defining, 


m 


x — x a 


g(s) = 


y~Va 


(3.3) 


x b -x a y b - y a 

generates x(s) and y(s) directly, once the function f and g are prescribed. 

Here we have used , 

f{s) = g{s ) = 5 (3.4) 


This procedure thus enables us in fixing the point distribution on the domain bound- 
aries. The interior grid points are obtained from a double interpolation as given by 
Fletcher [1] which we are not discussing here in this thesis. 
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One can see the point distribution on domain boundaries formulated bv 
stretching function as shown in Figure (3.3). 


3.3 Elliptic Grid generation 

* 3.3.1 Introduction 

An elliptic grid generation is a general method of generating boundary fitted coor- 
dinate system with Dirichlet boundary conditions on all boundaries. In case of an 
elliptic PDE grid generation technique, one coordinate is specified to be constant 
on each of the boundaries and a monotonic variation of the other coordinate around 
each boundaries is specified. Poisson’s equation in transformed plane of a linear 
elliptic system is used as a partial differential equation. Since only elliptic systems 
allow specification of a boundary conditions on the entirety of closed boundaries, it 
is used for any general closed boundary configurations. 


3.3.2 Co-ordinate System Generation 

Figure (3.1) shows the transformed plane of a grid for cascade domain. Transformed 
plane contains ^=constant and 77=constant lines. Hence our cascade domain is 
transformed from (x,y) to (f, tj) coordinate system. 

From the Figure (3.1), and T 2 are the solid boundaries, while T 3 and T 4 
are Inflow and Outflow boundaries. r 5 ,r 6 ,r 7 and F s are the periodic boundaries. 
The governing equations are elliptic partial differential equations and more specifi- 
cally the general Poisson’s equation is given by. 


+ £yy ~ > 'l) 


(3.5) 


CEN I r< !_ v v 


HASn 



3.3 Elliptic Grid generation 


52 


Tjxx ~b Viry — 


(3.6) 


where, P and Q are the known functions used to control interior grid clus- 
tering. 


i\ 


rs re T6 



$=o re 


n 


rs ^ =I 


% 


Figure 3.1: Transformed Plane of a domain 


With the Dirichlet boundary- conditions. 


£ 


& (x,y) 


m 


b(*,y) 


, [x, y] € Ti, 


, [x. y] 6 r 2 


£ 

V ) \ V 2 

and similarly for r 8 , T4 • • • and r 8 , w-here rji, 772 • • • and t]s are different constants 


{t)i < 772) while £2 1 • • and £ 8 are specified monotonic functions on Ti, T 2 • • • and 
r 8 respectively. 

For the numerical computation in uniform rectangular transformed plane 
all the variables have to be changed accordingly. Thus the transformed plane equa- 


tions are, 
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- 2 3xz„ + 7 ^ + J 2 (Pxz + Qx v ) = 0 
Q-Vtf ~ + ly-rm + J 2 (Piz + Qx v ) = 0 


where, 

a = x 2 + y 2 , 7 = X £ 2 + y ? 2 , 

0 X^X-q T y^y-r]i J — X^yr, y^Xrj. 

With the transformed boundary conditions, 



MtrVl) 


• [£> ^?i] £ Fi 


(3.7) 

(3-8) 



• [£> V 2 ] € Ta 


9 ifcrji) 

92 (^ 92 ) 

and similarly for T 3 , T 4 and r 8 , where the functions fi, f 2 ■ ■ ■ and g x . 
g 2 ■ ■ ■ are specified by the known shape of the contours Id, T 2 • • • T 8 and specific 
distribution of £ at these places. 

The equation (3.7) and equation (3.8) give the quasi-linear elliptic system 
for physical coordinate functions, x(£, 17 ) and y(£, rj) in a transformed plane. With 
a recommendation by Thompson [1], Poisson’s equation could be written as. 


£r* + £yy = ~ X a i S 9 n i^ - &)exp(-Ci | £ ~ fi I) 
1=1 
m 

- X bjSgn(£ - £j)exp{-dj((€ - £j) 2 + {rj- 9j) 2 ) 1/2 ) 
j - 1 

= p( «. 1 ) 

71 

Tj xx + g yy = - X a i s 9 n (9 ~ 9 ijexp(-Ci \rj-m\) 


(3.9) 
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- b J s 9n{r) - ^expi-djiiZ - ^j) 2 + (77 - Vj ) 2 ) 1 / 2 ) 


7 = 1 


Q( f,i?) 


(3.10) 


where the constants d£, &£, and d; are the grid lines controlling constants. In a 
given domain, close spacing between f = constant lines for any £ = & line or lines 
in a domain can be achieved by selecting first term from equation (3.9) as. 


P{£-.V) = - ^aisgn^ - £i)exp(-Ci | & |) ( 3 - 11 ) 

i=i 

The function sgn has the properties, 

sgn(x)=l if x is positive , 

=0 if x=0 , 

=-l if x is negative. 

Thus, the expression sgn(£ - &) is +1 for f > £ and -1 for f The 

constant a, is maximum magnitude of P and it controls the extent by which £- grid 
lines shift towards f The constant c, is a decay constant of P which basically 
controls the specific zone or region for which the grid to be attracted. The constants 
at and c t are not fixed for a physical system, but one has to get those values by trail 
and error until one gets the desired results. 

If instead of attracting towards a complete £=constant line, we w r ould like to 
attract £=constant line towards a given point (&, pi) in domain, the control function 
P could be reduced from equation (3.9) to, 

m 

P(£, 77 ) = - bjsgn{Z - ^j)exp(—dj((^ - £j) 2 + (j? - p j) 2 ) 1/2 ) (3.12) 

j = 1 

Similarly, choosing appropriate value of b{ and dj, we have to work out by' trial and 
error until we get satisfactory results. For the attraction to the line as well as to a 
point simulteniously, we have used full P and Q equations. 
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For attracting 77 -constant lins towards a given 77-constant line or a given 
point (f, 77 ) in domain, the similar treatment is to be followed in the equation( 3 . 10 ). 

A precaution is to be followed in the selection of the constant values of a, 
or bi that if they are too large, the Extremum principle may get violated and the 
grid lines may cross each other. For adaptive grid generation, the control functions 
can be selected according to the error of gradient of predicted solution. 

3.4 Results 

The grid generated with an algebraic PDE technique has been brought to desired 
level of spacing accuracy at the aerofoil solid boundary zone for the rj=constant 
lines with an elliptic PDE grid generation technique. By using equation(3.9), we 
get the spacing accuracy in 77 direction. As we have, aerofoils at 77=0 and 77=1 

in transformed plane, we are attracting the ? 7 =constant lines towards 77 , lines for 

* 

normal spacing of the order of 10 ~ 3 near the aerofoil surface. Figures (3.3) and (3.5) 
show the algebraically generated grids with spacing in 77 direction close to aerofoil 
is 10 ~ 2 , which is controlled by elliptically generating rhe same grid for the spacing 
accuracy in 77 direction as 10 -3 shown in Figures (3.4) and (3.6). 

The constants of a^s and c,-’ s are estimated by trial and error method which 
gives the desired level of spacing close to solid boundary of the order of 10 -3 . For 
the present code we have used constants as tabulated in Table (3.1). 


Grid 

a : 

bi 

cj | di 

193x179 

10.00 

15.00 

00.20 

00.09 

80x60 

05.00 

15.00 

00.10 

00.09 


Table 3.1: Values of the constants for different grids 
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\\ e have observed that for the cascade geometry, the constant values given 
above gives us the results shown in Figures (3.4) and (3.6). Any higher value of a* or 
d will cause violation of Extremum principle giving inter-crossing grid line domain. 

Due to sudden change in profile at leading and trailing edge, solution may 
diverge or there might be spurious oscillations of the solution. To avoid this, the 
smoothness of the grid lines near the leading edge and trailing edge is necessary. In 
the case of transonic or supersonic flows the shock wave generated at the leading edge 
while the expansion fan is found at trailing edge. To capture these shocks as well as 
expansion fans, we need to provide close grid spacing at leading and trailing edges as 
well as at £ iea( f=constant and £ tra j;=constant lines. In the code, we have forced the 
£=constant lines to be attracted towards £fead=constant and £trai/=constant lines by 
using equation (3.10). The point attraction is achieved by attracting 77=constant 
lines towards the (£,r])i ead and (^,T}) t raU- 









3.4 Results 


59 



Figure 3.4: Elliptic grid for NACA0012 straight cascade a. — 0° and 3 = 0° 
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'4 



Figure 3.5: Algebraic grid for DCA cascade a — 139.5° and 3 — 0° 
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Figure 3.6: Elliptic grid for DCA cascade a = 139.5° and 3 = 0° 


Chapter 4 


Results and Discussion 


A code for compressible flow is developed in FORTRAN to solve two dimen- 
sional unsteady viscous computations for flow past cascade. As a turbomachinerv 
application, the code has been developed for internal flow's which also can be used for 
any other physical domain for internal flows by choosing proper grid domain. The 
code is broadly divided into two parts. First, Euler code and second is the viscous 
code which adds the viscous terms along with a modifications in terms of boundary 
condition of the Euler code to give full Navier-Stokes solution. Since we have not 
used any local time stepping strategy, the results could be view r ed as time accurate 
inviscid and viscous solutions respectively. However, we have used an initial condi- 
tion which an uniform flow with u — and v = Wc everywhere. In that context 
the intermediate results may not be very meaningful. Also such an initial condi- 
tion will converge to actual steady state condition slowly. We recommend obtaining 
steady state Euler solution as a starting solution for Navier-Stokes computations. 
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4.1 Cases of Studies 

We have carried out our computations for the following cases, 

1. NACA0012 straight cascade: 

This case was considered as a test case for which some computational results 
are already available in [7]. 

In this cascade domain, we have used parallel aerofoils so as to have a 
straight computational domain, i.e. (a = 0° and 0 = 0°), where a is the inflow- 
angle and 3 is the outflow angle. The figure (3.2) shows the grid generated for 
straight cascade domain. All the computations for this configuration is carried out 
for a Mach number of 0.8. We have computed our solutions for two cases. First is 
for Euler computations with (80 x 60) grid domain and second is for Navier-Stokes 
computations for (193 x 179) grid domain. The grid spacing close to the solid wall 
for (80 x 60) domain is of the order of 10~ 2 , w-hile for (193 x 179) domain it is of 
the order of 10~ 3 . For Euler computations, v r e have used CFL number as 10 -3 . 
while Navier-Stokes computations are carried out for CFL number of 10 -5 . For 
Euler computations, we have results for the inflow-outflow 7 pressure ratio of 0.8. All 
Navier-Stokes computations are carried out for Reynolds number of 3xl0 5 w 7 hich is 
typical for compressor cascade operations. 

2. DCA cascade: 

This case of study is one of the cases documented in AGARD report [6]. 
For this case of study an experimental set-up is as follows: 


Aerofoil profile 

: DCA 

Number of blades (N) 

: 14 

Chord length (c) 

: 80 mm 

Aspect ratio (h/ c) 

: 3.75 
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Inflow angle (pi) 

139.5° 

Outflow angle (P 2 ) 

90° 

Blade turning angle (p s ) 

109.2° 


Xavier-Stokes equation is being solved for Reynolds number 5xl0 5 with a 
Mach number of 0.3. The grid has been generated for (193 x 179) domain with the 
solid wall spacing of 10 -3 . The computed solution is compared with the experimental 
results shown by Fottener [6]. The computations are carried out for CFL number of 
10 " 4 . 

4.1.1 Euler Solution for NACA0012 Cascade 

Euler computations are carried out to analyse the pressure variation along the aero- 
foil surfaces. Figures (4.1) to (4.4) give the pressure distribution calculated in terms 
of Coefficient of Pressure(C p ) along the upper and lower surfaces of the NACA0012 
aerofoil. These figures show the pressure distribution for various different time step 
levels. 

Figures (4.1) and (4.2) shows the C p variations for a inflow-outflow pressure 
ratio of 0.8. We could see that the solution after t = 1.0 has started capturing 
the weak shocks around third quarter of the chord. These shocks are due to the 
transients at that location. After t = 2.700, we could see that the sharp decrease in 
C p value from 0.6 to 0.2 around 0.65 of the cord length. It has been seen that the 
shocks are shifting from their previous locations as the iterations go on. The Figures 
at t = 4.5 to 5.0 give an idea about the shifting of this initial transient shock in 
downstre am direction. From the Figure at t = 10.0, we could see that the solution 
approached steady state as the shock move away from the cascade. The C p value 
variations from t = 12.2 to t = 19.8 show very little change, but as the solution 
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has not yet reached its convergence, we would see further changes in solution as 
iterations go on. It is to be emphasis that unlike the external flows, the internal 
flows approach their steady state value very slowly [8]. 

4.1.2 Navier-Stokes Solution for NACA0012 cascade 

For a CFL number of 10~ 5 . the pressure variation along the solid surface is shown 
in figures (4.3) to (4.8). It can be seen from these figures that, as the marching 
time step is very small, it takes longer time to show any significant changes between 
successive frames. 

From Figure (4.3). the results at t = 0.01 shows a primary oscillations 
around leading edge which are wiped out by t = 0.03 to produce smooth C p value 
variation. From t = 0.05 onward we could see a dip near the leading edge which 
strengthens as iterations proceed. The comparison of figures at t = 0.05 and t =. 0.15 
for C p shows significant change near the leading edge. The comparison of Figures 
between f = 0.15 and t = 0.25 shows that the dip has been shifted from 0.18 to 0.3 
of the cord length, while the maximum C v value is shifted from 0.98 to 0.95 of the 
cord length. We expect that the dip will continue to shift as well as steepen the C p 
value at 0.7 of the chord. 

4.1.3 Navier-Stokes Solution for DCA cascade 

For DCA cascade, Figures (4.9) to (4.17) show' the pressure variation in terms of C p 
along the aerofoil surface for M=0.3, w'hile Reynolds number is kept constant i.e. 
10 5 . 

From the Figures at t = 0.1 to t — 0.5, w r e could see that the initial transients 
in the solution show' rapid change in C p values. At t = 0.2, w r e could observe the 
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fluctuations at the leading edge for lower surface C p value, which retain its location 
up to t = 0.50. The upper surface C p , shows smooth variation from t = 0.1 to 
t = 0.2 which after t = 0.24. shows fluctuating up to t = 0.50. 

After t = 0.50, we can see that both upper and lower surface C p values have 
reached a maximum value at x/c=0.9. These maxima start shifting from trailing 
edge to leading edge as time progresses. The comparison of t = 0.61 and t = 0.95 
shows that the maximum C p values are shifted from 0.9 to 0.7 for upper surface and 
from 0.9 to 0.65 for lower surface. We could even observed that the maximum value 
of C p comes down as iterations proceed. Figures for t = 0.72 to t = 0.95 show that 
some wiggles are propagating from leading edge to trailing edge, which are pushing 
upper surface C p value to form a sharp peak as a weak shock. At t = 1.34. we could 
see the peak C p value which has shifted to half chord, while C p lower maximum has 
also shown similar shift to 0.3 of the chord. 

From the Figures t = 1.34 to t = 1.62 show that the C p lower values are 
smooth, while the shock on the upper surface moves toward leading edge. The com- 
parison between t = 1.66 and t = 1.98 shows that the C p upper value curve displays 
weakening of the shock while shifting towards the leading edge. The maximum C p 
upper value is reducing as the iterations proceed. We could also see that the C p 
lower value curve is also smooth with the minimum value shifted from 0.4 of the 
chord to 0.2 of the chord. After every 100 iterations we could see the comparable 
changes in the solution. 

As this computations are not yet over, we could expect from the solution 
pattern that the C p lower value will go further down to get minimum value at 0.7 of 
the chord, while C p upper will show the maximum value at around 0.2 of the chord 
with the w r eak shock pattern. 
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4.2 Conclusions 

1. Grids for the cases presented are generated by an algebraic as well as 
an elliptic grid generation technique. Although, an algebraic grid is much faster 
to generate, an elliptic grid gives smoother variations for J and hence a faster 
convergence. To capture accurate solutions we should have very fine grid spacing 
near the solid boundary. For our cases of study, we should use much more refined 
grid close to the solid surface in comparison with the present grid spacing. 

2. We have used higher order interpolation upwinding scheme in calculating 
flux matrices. The spurious variations in the solution for higher Mach numbers could 
be minimised by using nonlinear Flux-Limitter(FL) approach [20]. This flux-limitter 
controls the Gibb's phenomenon generated for higher wave numbers^). 

3. The computations are carried out for straight NACA0012 cascade as well 
as DCA cascade. From the results we could see that the straight cascade is more 
difficult to compute compared to DCA cascade, since the variations in the solution 
are rather slow. For the DCA cascade study case the results are compared with the 
experimental [6] case study. As the computations are still going on, we could expect 
from the results so far produced, the solutions should approach to the experimental 
results. 


4.3 Future Scope of the Work 

We have studied and analysed the results for straight as well as angular cascade 
domain with transonic as well as subsonic inflow Mach numbers. 

In future, we would like to investigate solution for the high Mach number 
inflows at transonic or supersonic speeds. In the present work, our outflow boundary 
conditions are predefined, but in actual flows, the outflow' velocity angles define the 
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runiet crs at outflow. This feature can be incorporated in our code 

"“CVeuact simulation of the turbine flow. We would like to study following 
to actu^ e * 

K ts regarding the cascade flows in future: 

, Noise added to inflow to study effect of free-stream turbulence. 

2 Tate the outflow solution of the present set of computatron as an mflow 

,, pr computation to investigate effect of gust. • _ 

“ 3. Here periodicity is over row of cascades- it can be extended for tnuluple 

. ctiidv events like rotating stall. 

"" tee we are obtaining time accurate solution, the present set of resnlts 

.betlL ^yrng — * - — aspect of cascade flows. 
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Figure 4.3: Pressure Variation for NACA0012 Navier-Stokes Solution Cont n . Mach 
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Figure 4.9: Pressure Variation for NACA0012 Navier- 
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Figure 4.14: Pressure Variation for DC A Navier-Stokes Solution Cont 
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